On the adaptive design rules of biochemical networks in evolution.

Biochemical networks are the backbones of physiological systems of organisms. Therefore, a biochemical network should be sufficiently robust (not sensitive) to tolerate genetic mutations and environmental changes in the evolutionary process. In this study, based on the robustness and sensitivity criteria of biochemical networks, the adaptive design rules are developed for natural selection in the evolutionary process. This will provide insights into the robust adaptive mechanism of biochemical networks in the evolutionary process.We find that if a mutated biochemical network satisfies the robustness and sensitivity criteria of natural selection, there is a high probability for the biochemical network to prevail during natural selection in the evolutionary process. Since there are various mutated biochemical networks that can satisfy these criteria but have some differences in phenotype, the biochemical networks increase their diversities in the evolutionary process. The robustness of a biochemical network enables co-option so that new phenotypes can be generated in evolution. The proposed robust adaptive design rules of natural selection gain much insight into the evolutionary mechanism and provide a systematic robust biochemical circuit design method of biochemical networks for biotechnological and therapeutic purposes in the future.


Introduction
Robustness is a ubiquitously observed property of biological systems. It is considered to be a fundamental feature of complex evolvable systems. It is pointed out that robustness facilitates evolvability and robust traits are often selected by evolution (Kitano, 2004), i.e. complex biological systems must be robust against environmental and genetic perturbations to be evolvable. Evolution often selects traits that might enhance robustness of the organism.
The central role of biochemical networks in cellular function provides a strong motivation to search for the underlying principles of adaptive evolution of biochemical networks. In this study, in order to test whether a physiological function would prevail under a new environment or not, the robustness and sensitivity criteria are developed to measure the tolerance of the metabolite concentration values of a biochemical network in the face of environmental changes. That is, we derive necessary and sufficient conditions for the metabolite network to be preserved by natural selection in the evolutionary process.
The evolutionary analysis is based on two concepts, natural selection and evolution (Freeman and Herron, 2001). In the past, most molecular biologists and biochemists assumed that variations in biochemical networks were mainly due to historical accidents and natural selection. But the design principles of biochemical networks via natural selection in evolution are still in conceptual description, not in mathematical rules. Can these mathematical natural selection principles for biochemical networks in evolution be unraveled? The investigation of design principles of biochemical networks is in its infancy and more underlying rules remain to be discovered. In fact, robustness allows change in the structure and components of the system owing to these perturbations and disturbances, but specifi c functions are maintained. Hence, robustness facilitates evolvability and evolution selects robust traits (Yi et al. 2000;Kitano, 2004).
In this study, the robustness criterion based on S-system is derived as a necessary adaptive design rule of a biochemical network under natural selection. In addition, to guarantee small changes of metabolite concentration values under environmental disturbances in the evolutionary process, some sensitivity criteria are proposed as the sufficient conditions for the adaptation of a biochemical network so that it can play its proper role in the corresponding physiological system. Two adaptive design schemes of robustness improvement are developed for biochemical network evolution. One scheme is to compensate the effect of parameter variations to meet the robustness criterion easily. In this way, some redundant and selfregulatory pathways are selected by natural selection to attenuate the effect of parameter variations. The other is to enhance the system structure stability of biochemical networks to tolerate larger parameter perturbations. In this way, negative feedbacks and positive feedbacks are selected to improve structure stability. These two adaptive schemes are two design methods of biochemical networks to improve their robustness and to maintain the live function against environmental changes in the evolutionary process. The biochemical networks with improved robustness can survive under natural selection. At the same time, the sensitivities of biochemical networks to environmental disturbances are also attenuated to maintain their metabolite mechanisms for normal physiology, i.e. insusceptible to environmental disturbances for species evolution or diseases such as virus infection, with immunity for individual. They are also considered in the adaptive design rules of natural selection in evolution.
Since many solutions can meet the robustness and sensitivity criteria of natural selection, a variety of biochemical networks may survive in evolution. A variety of biochemical networks with some structural differences may arise in evolution. For example, in the TCA cycles in different species, their fi nal products are almost the same from yeast to human, but their biochemical networks have some structural differences in intermediary biochemistry reactions. Based on the adaptive design rules of biochemical networks via natural selection, one possibility of the diversity in biochemical networks in the evolutionary process is to increase the complexity of networks through successive addition of feedback and feedforward pathways to enhance robustness against genetic mutations and environmental perturbations (Barkai and Leibler, 1997;Alon et al. 1999;West-Eberhard, 2003).
Natural selection can select only from the mutated biochemical networks that already exist in nature and cannot instantly create a new and optimal biochemical network (or phenotype) to maintain the live function. The co-option of existing biochemical networks to new networks is one of the crucial features in evolution (Kitano, 2004). Several biochemical networks are combined through positive feedback loops and negative feedback loops so that normal cellular physiology and developmental processes can be maintained. This intrinsic robustness of a biochemical network enables co-option, so that new morphologies can be generated in the evolutionary process (Freeman, 2000;Kitano, 2004). The proposed adaptive design rule via natural selection can mimic the evolution of biochemical networks by computational simulation. Simply speaking, the evolutionary way is to improve its robustness of biochemical networks to tolerate the parameter variations and environmental variations to play their proper role in the corresponding physiological systems. The proposed robust adaptive design rules of natural selection also provide a systematic robust biochemical circuit design method of biochemical networks for drug design and robust engineered synthetic biocircuit design purposes in the future. We can use computational prediction and rational design (Altamirano et al. 2000;Johannes et al. 2005;Tsuji et al. 2006), directed evolution (May et al. 2000;Wang et al. 2000;Yuan et al. 2005) and dynamic controller (Farmer and Liao, 2000;Bulter et al. 2004) to quickly create a library of variants for artificial evolution to achieve the desired property of biochemical networks. Rational design and directed evolution are to modify the catalytic or binding property of an enzyme which corresponds to the changes of kinetic parameters g ij and h ij in the S-System model through modulating the enzyme structure and through DNA shuffling, respectively. A dynamic controller is to construct a feedback or feedforword pathway. Finally, a computational simulation example is given to illustrate the adaptive design mechanism of a biochemical network via natural selection in evolution.

Notations
For a vector x = [x 1 , …, x n ], the l 2 norm of x is defi ned as ||x|| 2 = x x x n 1 2 2 2 2 + + + . We say x ∈ l 2 , if ||x|| 2 < ∞. For a matrix A and y = Ax, the l 2induced matrix norm is defi ned as ||A|| 2 = sup x∈l2 y x 2 2 , i.e. the gain from x to y. It has been shown that ||A|| 2 = σ max (A) = max i λ i T A A ( ), where σ max (A) denotes the largest singular value of A and λ i (A T A) denotes the ith eigenvalue of A T A. ||A|| 2 < 1 if and only if AA T < I i.e. A is contractive, where I is the identity matrix (Gill et al. 1991;Weinmann, 1991).
The mathematical model, robustness and sensitivity analyses of a biochemical network under natural selection in the evolutionary process are introduced at fi rst.

Model of a Biochemical Network
In a biochemical network, one often measures rates of reaction or infl ux and outfl ux rates of substrates, enzymes, factors or products and the rates correspond directly to changes in concentrations. The S-system model has been developed to write the reaction relationship of metabolites in differential equations in terms of their concentrations. The dynamic system of a biochemical network is described in the following S-system representation (Voit, 2000;Savageau, 2001)  where X 1 , …, X n+m are metabolites, such as substrates, enzymes, factors or products of a biochemical network, in which X 1 , …, X n denote n dependent variables (intermediate metabolites and products) and X n+1 , …, X n+m denote the independent variables (initial reactants and enzymes). i X , the rate of change in X i , represents concentration change of a dependent variable due to production (accumulation) or degradation (clearance). Each term is the product of the rate constant, α i or β i , which is positive or zero and all dependent and independent variables that affect directly the production and degradation reaction, respectively.
Each variable X j is raised to the power of a kinetic parameter g ij and h ij , which represents that X j activates (inhibits) X i when its value is positive (negative). The rate constants α i and β i and kinetic parameters g ij and h ij are related to the characteristics of the biochemical network. The nonlinear Equation (1) describes the dynamic evolution among dependent variables. How to construct the S-system representation of a biochemical network and how to estimate its parameters from experimental data can be found in the classic textbooks (Savageau, 1976;Voit, 2000) and references therein. Recently, the nonlinear parameter estimation problem of Ssystems has been efficiently solved by evolution optimization methods (Tsai and Wang, 2005;Ko et al. 2006).
Measuring directly the robustness of the nonlinear system in Equation (1) is difficult most of the time. Fortunately, the phenotype (metabolite concentration values) of a biochemical network is close to the steady state, i.e. the transient time to dynamic equilibrium is short enough in the real world and the steady state of biochemical networks can be analyzed by simple algebraic methods. Therefore, we shall focus on the robustness of a biochemical network at steady state in this paper.
Consider the steady state of biochemical network in Equation (1), i.e. the production and degradation of each dependent variable is balanced (Voit, 2000). (1) After some rearrangements, we get Introduce new variables and coefficients as follows: The steady state of a biochemical system is written in n linear equations in terms of n + m variables as follows (Voit, 2000) a In Equation (6), the dependent variables are separated from the independent variables. Let us denote where A D denotes the system matrix of the catalytic interactions among dependent variables and A 1 indicates the catalytic interactions between the dependent variables Y D and the independent variables Y 1 (i.e. environmental medium to the metabolic system). We then obtain the steady-state equation in the nominal parameter case as follows

Natural Selection Criteria for Biochemical Networks in Evolution
From Equation (7), if the inverse of A D exists, then Y D can be solved uniquely. It means that the biochemical network will result in only one steady state as long as A D −1 exists. The assumption makes sense and agrees with the real biochemical networks. The steady state (or phenotype) of the biochemical network is solved as follows (Voit, 2000): Biochemical systems perform their physiological function within some local region in system parameter space in the evolutionary process. They tend to be robust to local changes in the values of the parameters that defi ne the system in Equation (8). In the evolutionary process, suppose that some parameter variations Δα, Δβ, Δh, Δg and ΔY I , which could be considered as design parameters in the evolutionary process owing to genetic mutations or environmental changes, alter the kinetic properties of a biochemical network in comparison with the nominal kinetic parameter case in Equation (7) as follows: where the parameter variations of the biochemical network are defi ned by   The robustness analysis of a biochemical network in this study is to check the tolerance for kinetic parameter variations with respect to the maintenance of normal physiological function of perturbed biochemical networks in the evolutionary process. First, the ΔA D will infl uence the existence of the steady state. Then, the variations Δb, ΔA I and ΔY I will infl uence the sensitivity of biochemical networks to the environmental variations in the evolutionary process. The effects of these parameter variations (i.e. the design parameter space in evolution) on the biochemical network at the steady state (i.e. i X = 0) will be discussed in the following paragraphs.
Equation (9) is equivalent to By the similar analysis from Equation (7) to Equation (8), one can show that the condition that the system in Equation (10) can be solved uniquely is the existence of the inverse of (I + A D −1 ΔA D ). It has been shown that if the following robustness criterion holds (Gill et al. 1991;Weinmann, 1991;Nobel and Daniel, 1998;Chen et al. 2005) 1 exists and the phenotype (steady state) of the perturbed biochemical network in Equation (10) is uniquely solved as follows The physical meaning of Equation (11) and (12) is that if the inverse ( ) 1 exists, the phenotype can be preserved with some variation under this parameter variation ΔA D , i.e. if the robustness criterion in Equation (11) is satisfi ed, the parameter variations ΔA D can be tolerated by the system structure of the biochemical network A D in the evolutionary process and the biochemical network tends to be robust to local changes in the values of parameters that defi ne the system. Otherwise, the parameters reach a threshold beyond which the behavior of the biochemical network changes dramatically and the phenotype may cease to exist, i.e. the individuals with parameter variations (design parameters in evolution) that violate the robustness criterion in Equation (11) will be eliminated by natural selection. Therefore, the perturbed biochemical network should satisfy the robustness criterion in order to guarantee the existence of its dynamic equilibrium (for the normal physiological function) in the evolutionary process. Because the violation of Equation (11) means a lethal perturbation, it is the necessary condition to survive under natural selection. From the robustness criterion in Equation (11) are so much redundancy due to duplicated genes, modularity, self-regulation and feedback circuits in the biochemical networks in nature (Isaacs et al. 2003;Langkjaer et al. 2003;Kellis et al. 2004;Teichman and Babu, 2004).
However, the satisfaction of the robustness criterion in Equation (11), i.e. the parameter variations T is bounded by the system structure matrix A A D D T , does not always mean the perturbed biochemical network will survive in evolution because it only guarantees the existence of the steady state. But the phenotype (steady state) may be far from the nominal value for the normal physiological function. In order to play its proper role in the corresponding physiological system, its metabolite concentration values should not change too much from the nominal value. In this situation, the biochemical network should be less sensitive to the other parameter variations and environmental changes. This is the suffi cient condition for natural selection for a perturbed biochemical network to survive under natural selection. In the above robust analysis, we only discussed the effect of kinetic parameter variations ΔA D on the existence of the steady state of a biochemical network. Now, let us consider the sensitivities to the variations of the other parameters Δb, ΔA I and the change of the environment ΔY I in the evolutionary process.
The changes Δb, ΔY I , and ΔA I will infl uence the variations of steady states Y D . Their effects on Y D have been discussed by the sensitivity analysis of biochemical network (Savageau, 1971;Ni and Savageau, 1996ab;Voit, 2000), i.e.
In order to tolerate the variations Δb, ΔY I and ΔA I to preserve the phenotype of the biochemical network in the evolutionary process, the sensitivities in Equation (13) should be below some values as follows where s 1 , s 2 and s 3 are some small sensitivity values so that the phenotypes of perturbed biochemical networks would not change too much in comparison with the nominal values in Equation (13) and can be favored by natural selection, i.e. the sensitivity criterion in Equation (14) can be considered as the suffi cient condition of natural selection for biochemical network evolution. In general, the sensitivities s 1 , s 2 and s 3 are chosen as the sensitivities at the nominal case, because the nominal (healthy) biochemical network is less-sensitive to parameter variations and environmental changes. Based on Equation (13), Equation (14) can be written in the following equivalent form,

I s A A A A s A A Y Y s A
That is, Equation (15) determines the ranges of the sensitivities of Y D to parameter variations and environmental changes by natural selection in the evolutionary process. For a functional biochemical network, it should satisfy the sensitivity criteria in Equation (15) to confi ne the metabolite concentration values not to be changed too much. Hence, the steady state (phenotype) of a biochemical network can be preserved while exposing the parameter variations and environmental changes to natural selection in the evolutionary process. This can be considered as a suffi cient condition for survival for the biochemical network.
The robustness criterion in Equation (11) and the sensitivity criterion in Equation (15) are together considered as the criteria of natural selection in evolution. If one of them is violated, it will lead to the dysfunction of the biochemical network and the perturbed biochemical network will be eliminated by natural selection. Therefore, the robustness criterion in Equation (11) and the sensitivity criterion in Equation (15) could be considered as the adaptive design rules of biochemical networks by natural selection in the evolutionary process. The specifi cations of sensitivities s i , i = 1, 2, 3 in Equation (15) are species by species. In general, these sensitivities should be small in order to avoid too much infl uence from the environmental disturbances in the evolutionary process. ) and the steady state will cease to exist. However, the equality could hold in Equation (14) because we do not want the sensitivities of perturbed systems to be larger than the sensitivities of the nominal system, which has no singular problem. (ii)Actually, the sensitivity matrices in Equation (13) hold if all the perturbations Δb, ΔY I , ΔA I are very small (Voit, 2000) (it was originally derived by . For the convenience of discussion on perturbations, it was modifi ed to the form in Equation (13)). If some perturbations are large, the equalities may be violated. One proposition of Theory of Evolution is that "Gradual evolution results from small genetic changes that are acted upon by natural selection" (Freeman and Herron, 2001). Obviously, in evolutionary process, Δb, ΔY I , ΔA I are all assumed to be small in every change. In this situation, the equalities in Equation (13) always hold. (iii)The assumption that the three sensitivity inequalities in Equation (14) all hold for natural selection is based on the fact that biochemical networks are the backbones of physiological systems and can not be too sensitive to environmental changes especially for some core (conserved) biochemical networks. If some sensitivity criteria in Equation (14) are relaxed, i.e. some of inequalities in Equation (14) are violated, the phenotypes with changes to some environmental variation will also be favored by natural selection. In this situation, the phenotypes of biochemical networks are much infl uenced by environmental variation that they may be more adaptive to the environmental changes in the evolutionary process. In this case, new phenotypes are more easily generated to adapt the new environment. They will be discussed in the sequel.

Computational Examples
An example is given below to illustrate the mathematical adaptive design rules of biochemical networks by natural selection in the evolutionary process. Consider the following biochemical network (Savageau, 1976;Voit, 2000).
The biochemical network in Equation (16)  Y I || 2 . That is to say, the perturbed networks to be selected by natural selection should have less sensitivities than the nominal biochemical network. By the adaptive design rules based on robustness and sensitivity criteria, the perturbed biochemical network (I) in Equation (19) violates the robustness criterion in Equation (11) and the parameter variations ΔA D1 can not be tolerated by the biochemical network. In this situation, the biochemical network (I) will be eliminated by natural selection without consideration of sensitivities. More precisely, the set of parameter variations ΔA D 1 due to mutations is lethal. Though the biochemical network (II) in Equation (20) satisfi es the robustness criterion in Equation (11), its steady state is farther from nominal value and violates sensitivity criterion in Equation (15). It means that the phenotype of biochemical network (II) is easier (more sensitive) to be destroyed while exposing to environmental disturbances, i.e. ΔA I , Δb or ΔY I due to environmental changes. Hence, there is a large probability that the biochemical network (II) will be eliminated by natural selection. More precisely, the parameter variations due to mutations ΔA D2 are not lethal, but the biochemical network is susceptible to environmental disturbances for species evolution or diseases such as virus infection, i.e. less immunity for individuals.
From Figure 2, it is seen that the steady states, Y D3 + ΔY D3 and Y D4 + ΔY D4 , of biochemical network (III) in Equation (21) and biochemical network (IV) in Equation (22), respectively, are all close to the nominal values of the steady state in the nominal biochemical network of Equation (16) in Figure 1. In addition, the robustness criterion in Equation (11) and the sensitivity criterion in Equation (15) are all satisfi ed so that the variations ΔA D3 and ΔA D4 , Δb 3 and Δb 4 , ΔA I3 and ΔA I4 as well as environmental disturbances ΔY I3 and ΔY I4 do not affect the normal function of the biochemical networks too much. In other words, the biochemical networks (III) and (IV) are robust to intrinsic parameter variations and less sensitive to environmental variations, so that the two biochemical networks must be more favored by natural selection. In the next generation, the other perturbed biochemical networks will be selected by natural selection with the same procedure. This co-option of existing biochemical networks to new networks by natural selection is considered one of the crucial features in the evolutionary procedure. Several biochemical networks combined with negative and positive feedback loops are robust against parameter variations and environmental disturbances so that normal cellular physiological and developmental processes can be maintained. This intrinsic robustness and sensitivity of the biochemical network enable co-option, so that the new phenotypes can be generated by natural selection (Kitano, 2004). Therefore, the robustness criterion in Equation (11) and sensitivity criterion in Equation (15) can be viewed as the mathematical adaptive design rules of natural selection.
From Table 1, we can fi nd that the sensitivities of biochemical network (III) and biochemical network (IV) are both smaller than the nominal one. That is, there is a high probability that the mutated biochemical networks with smaller sensitivities s 1 ,  Table 1.

Diversity of Biochemical Network within Organisms or Individuals in Evolution
There are many perturbed biochemical networks that can satisfy the adaptive design rules of natural selection in Equations (11) and (15) to survive in the evolutionary process. If they are all selected by natural selection, there will be some differences in phenotype (see Equation (12)) among these selected biochemical networks. However, as the values of parameters continue to change, they reach a threshold (i.e. the robustness criterion in Equation (11) is violated) beyond which the behavior of the biochemical network changes dramatically. It may thus settle in a new local region of another steady state with a different set of behaviors, or it may become completely dysfunctional and not survive under natural selection in evolution.
After several generations in the evolutionary process, due to co-option of existing biochemical networks to new networks, diversities of biochemical networks with conserved physiological function but with different structures will be developed (Freeman, 2000). However, if the requirements on the robustness in Equations (11) and the sensitivity in Equations (15) are more strict (or more conservative), only few solutions (or structures) can be selected by the natural selection to meet these requirements. This is the reason why a conserved core biochemical network has less  (25) and its time responses. The biochemical network V will adapt to the environment with large variation ΔY I in evolution. The biochemical network will adapt to an environment with large variation in Δb and ΔA I . diversity (Kitano, 2004). For examples, in the evolutionary process of the TCA cycle, the pentosephosphate pathway and the glycolysis pathway in different species, their fi nal products are almost the same from yeast to human. However, their biochemical networks have some differences in intermediary biochemistry reactions. From the simulation examples of perturbed biochemical networks in Equation (19), (20), (21) and (22), the perturbed biochemical networks (III) and (IV) in Equations (21) and (22), which are shown in Figure 2c and 2d, respectively, can be seen as the diversities of the biochemical network in the evolutionary process.

Remark 2
Since the violation of sensitivity criteria of natural selection in Equation (15) is not lethal, the relaxation of some sensitivity inequalities in natural selection will make biochemical networks more easily adapt to new environmental changes. For example in Figure 3a and 3b, suppose the biochemical network suffers ΔA D5 and ΔA D6 , respectively, due to genetic mutations in the evolutionary process in Equation (23), the biochemical network in Equation (16) is then perturbed to the (V) and (VI) network in Equation (24) and (25), respectively. In the former case, the second sensitivity criterion in Equation (15) is violated. The biochemical network will adapt to the environment with large variation ΔY I in evolution. In the latter case, the fi rst and third sensitivity criteria are relaxed.
The biochemical network will adapt to an environment with large variation in Δb and ΔA I . The sensitivities of biochemical networks in Equation (24) and (25) are listed in Table 2.

Conclusion
From the biochemical network evolution point of view, a biochemical network should be robust enough to maintain its proper role in the physiological system under parameter variations due to mutations and environmental disturbances. On the basis of stability robustness and less sensitivity to the effects of genetic mutations and environmental variations, the proposed design rules (robustness and sensitivity criteria) are developed as the underlying mathematical adaptive design principles of natural selection in the biochemical network evolution. That is, in the evolutionary process, organism enhances the structure stability by feedback or feedforward circuits to improve the robustness of a biochemical network to tolerate parameter variations or compensates the effect of external or internal perturbations. The self-regulation and redundancy are of this kind of robust design favored by natural selection in the evolutionary process (Becskei and Serrano, 2000;Issacs et al. 2003). Therefore, in the evolutionary process of biochemical networks, robustness is the maintenance of specifi c functionalities of the network against perturbations, and it often requires the biochemical network to change its mode of operation in a fl exible way. In other words, robustness allows changes in the structure and components of Biochemical network in Equation (16) 3.4191 2.6647 2.5643 Biochemical network in Figure 3a 3.0248 3 2.2686 Biochemical network in Figure 3b 4.1534 1.4967 3.1151